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Abstract 

It has been suggested in the literature that it may be possible to locate the QCD critical end 
point using the Taylor series of thermodynamic variables about the u = axis. Since the phase 
transition at the critical end point is believed to be in the 3D Ising universality class, it would 
seem natural to test this method with the Ising Model, for which the answer is already known. 
The finding is that it is in fact possible to pinpoint the location of the Ising critical point using 
Taylor coefficients. 



I. INTRODUCTION 



The existence and location of a critical end point (CEP) on the QCD phase diagram is 
an issue of active debate in the theoretical [IHE] and experimental [9TTT2"] physics community. 
Perturbation theory tells us that at sufficiently high temperatures QCD matter exists in a 
deconfined quark gluon plasma (QGP) phase [T3HT5] . Since we are well familiar with the 
low temperature hadronic matter phase [16J, at some point on the /i — T plane a phase 
transition or at least a crossover must occur. 



The story is as follows, which begins with Fig. 1(a) a plausible sketch of the QCD phase 
diagram for physical values of strange, up and down quark masses. The generic features are 
a first order phase boundary that terminates at a CEP (/i c , T c ), and a crossover between the 
QGP and hadronic phases at \i = 0. A brief remark about the notation: T c refers to the 
temperature at the CEP, while T c / Q denotes the QCD crossover temperature at fi = 0. 

A first order phase boundary that terminates at a CEP is known to exist for exactly 
solvable random matrix models that have the same symmetries as QCD [T7j. Furthermore, 
the data that we have obtained from lattice QCD [HI [19] suggests that we should expect to 
see a crossover to the QGP along the \i = axis at T c / Q ~ 150 MeV [20]. Unfortunately, it 
is very difficult to test these predictions experimentally, hence determining the exact layout 
of the QCD phase diagram remains an open problem. 



For instance, in place of Fig. 1(a) , we may have a phase diagram with a first order curve 
that continues all the way to the fi = axis. Or, alternatively, a first order phase boundary 
and CEP may simply not even exist. It is known theoretically [21] that the transition 
(whether it be first order, second order or simply a crossover) to the QGP at fi = is 
heavily dependent upon, amongst other things, the strange quark mass m s . 

Assuming that the QCD phase diagram is that as depicted in Fig. |l(a)[ then it would 
be of great interest to make a theoretical prediction about the location of the CEP. In 
principle, this can be achieved via Lattice QCD by a Monte Carlo sampling of the QCD 
partition function 

Zqcd(T, {iM f }) = I VUe-^HVetMfifif), (1.1) 
J f 

however, it is well known that at non-zero chemical potential one encounters the notorious 
fermion sign problem. 

In spite of this limitation, the location of the CEP may still be potentially determined 
by lattice methods [22H21] if one performs a Taylor expansion of thermodynamic variables 
about fi = (or more generally, Re(/x) = 0). Since all of the Taylor coefficients are evaluated 
at zero or purely imaginary fi, there is no fermion sign problem. 

Since thermodynamic variables (pressure and its derivatives for instance) are non-analytic 
at the CEP, the radius of convergence of such a Taylor series specifies the exact distance 
from the series' origin to the nearest non- analytic point. If this point happens to be on the 
real [i axis, one has effectively located the CEP. 

The main drawback is that in practice, one is only able to compute Taylor coefficients 
up to a specific finite order. Locating the CEP, nevertheless, requires knowledge of the full 
Taylor series. One would hope that it is sufficiently reliable to simply extrapolate from the 
few knowable lowest order coefficients. 

The focus of this work is to test the effectiveness of the Taylor series method. If it has any 
hope of predicting the location of the CEP in QCD, it should at least be able to do so for 
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a much simpler toy model 1 whose second order phase transition is in the same universality 
class. We are, of course, talking about the 3D Ising Model, whose phase diagram is shown 
in Fig. |T(bjj 



The 3D Ising Model is not exactly solvable, but numerical simulations are able to yield 
high precision estimates of its critical exponents and critical temperature T c (unlike for 
QCD). I.e., the location of the Ising Model critical point is known. 

Looking at Fig. [2j the series expansion method in QCD amounts to walking along the 
fi = axis (where there is no sign problem) and computing the Taylor coefficients of the 
quark number susceptibility (QNS) at more or less evenly spaced intervals of T. It is from 
these coefficients that one tries to learn about what takes place at non-zero fi. The direct 



analogy, in the Ising Model, is to walk along some curve on Fig. |l(b) , and compute a similar 



set of Taylor coefficients (again at roughly even intervals) . However, since there is no sign 



problem anywhere on Fig. 1(b), we are free to pick a convenient shape for this curve. For 
the sake of simplicity, we will consider the problem of locating the critical point based on 
correlation functions measured along a line parallel to the H = axis, as depicted in Fig. [2} 

Having computed the Taylor coefficients of the Ising Model free energy F up to some given 
order, we will then analyze the properties of this truncated Taylor series in an attempt to 
locate the critical point. Specifically, since the magnetic susceptibility diverges at the critical 
point, we are looking for a non-analytic point of F. For the Ising Model, this is a very 
convoluted way of locating the critical point: but that is the point. We want to test how 
the extrapolation method works in an example where we can also directly locate the critical 
endpoint, to check our work. 

The point is to mimic what has been done or at least is presently doable for QCD. Due 
to the ease with which one can perform statistical simulations of the Ising Model we are 
theoretically able to compute the Taylor coefficients needed for this analysis up to some very 
high order, so we are not restricted in this regard. In QCD, this is definitely not the case. 
Hence, we really only want to go to high enough order to at least address the issue of what 
is needed to learn anything. 



II. TAYLOR SERIES ABOUT H = H VS. TAYLOR SERIES ABOUT fi = 

For an overview of the Taylor series method applied to QCD the reader is referred to 
|24j . Furthermore, a detailed overview of the Ising Model is given in [26]. We will now show 
how this translates over to the Ising Model, which, in cartoon form, is illustrated in Fig. [2j 
Assuming that the QCD free energy is extensive, the thermodynamic pressure is given by 

T 

Pqcd = ylogZ QCD (2.1) 

which can be expanded in a dual Taylor series in the up/down quark chemical potentials 
about fi u = = 

P QCD (T^ u ^ d ) = P QCD (T,0,0) + Xn^n d —^. (2.2) 

riu,n d 



Svc also 
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FIG. 1: On the left, an abridged version of the QCD phase diagram for physically relevant values 
of m s and m UjC i, showing the existence of a CEP. On the right, Ising Model phase diagram when 
D > 2. 



This series can be reformulated in terms of baryon = + /x^)) and isovector (/ij = 
\{Hu — fJ'd)) chemical potentials. Fixing \x u = fid and denning fi = fib/ 3 leaves a simplified 
Taylor series in fi 2 = = fi d . 

The various Xi,/ S have expressions in terms of the expectation values of products of My, 
Mj l and its derivatives. Moreover, since the origin of the Taylor series is fi = 0, all of these 
computations are performed with DetMy real and non-negative. 

To construct the 3D Ising Model analogue, we begin by denning an N x N x N lattice E 
as a set of individual spins <Ji on a 3D simple cubic lattice with volume V = N 3 and lattice 
spacing a = 1 

V = {*!,.. .,<r v }. (2.3) 
These spins interact via the Hamiltonian 

i <i,i> 

where the second summation is over all nearest neighbor pairs, and J > is ferromagnetic. 
Given H, an expression for the partition function then follows. Defining the total spin a 

ff=^X>* (2.5) 

we have 

Z = tr{e"^} = ]T eP HV ° J] (2.6) 

s (ij> 

where, as usual, (3 = 1/T. Hence, the free energy per unit volume is 

F = -^-logZ (2.7) 
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(a) Ising (b) QCD 

FIG. 2: Visualization of the series expansion method described in |24j and this work. Circles of 
radius r ro about any given expansion point are depicted in green. Domains of convergence may 
not necessarily be finite sized for general T ^ T c (i.e. r^T ^ T c ) = oo); however, when r^T) 
is finite, the corresponding singular point must be complex. Except, of course, when T = T c , at 
which point the singular point must be strictly real. 



from which we obtain the thermodynamic variables 

- - -(ml w> 

* - -Or M 

where M and \ are respectively the magnetization and magnetic susceptibility. A Taylor 
series for F is obtained by simply noting that Z is a moment generating function, and hence 
log^ is a cumulant generating function [27] . Expressions for the various cumulants K n in 
terms of the moments (a n ) can be obtained from the generating function 

n n = —\og(e^) (2.10) 



dt n 



t=Q 



and, for instance, 



K X = (a) 

k 2 = (a 2 ) - (a) 2 

k 3 = (a 3 )-3(a 2 )(a) + 2(a) 3 

k 4 = ... (2.11) 
Thus, the Taylor series for F is given by 

F(H;H ,T) = - v log Z(H ,T) - ^E^vVv (2 " 12) 
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which is an expansion in H at fixed T. 

In summary, the analysis for QCD and the Ising Model amounts to extrapolating the 
behavior of the higher order Xi,/ S an d the /%'s from those that can be calculated. If a 
singularity is suspect, its distance from the expansion point is equal to the series' radius of 
convergence. For an arbitrary function f(x) expanded in a series 



/(*) = !> 



X ' 



(2.13) 



n=0 



the radius of convergence r^, is given formally by 



limSUp,,^ VlTn 

In our case, we may infer the radius of convergence from the asymptotic behaviour of 

fn 



(2.14) 



ra+1 



Alternatively, we may also consider the n — > oo limit of 



2n 



2n+2 



or 



Pn 



/o 



in 



(2.15) 



(2.16) 



(2.17) 



All of these definitions are effectively equivalent in that they all yield the same n — > oo 
limit. However, r n and p n are of relevance to [21] since the Taylor series that they analyze 
only has even powers of /i. This is, however, a minor technicality. In both models, the 
asymptotic behavior of r n (or r n or p n ) is inferred from the series' few lowest order Taylor 
coefficients. Specifically, the highest order that is achieved for QCD is n = 8. So that we 
remain on roughly the same footing, we will thus restrict our analysis to n < 8 coefficients. 

Monte Carlo studies [28] with lattices up to iV = 256 have revealed the Ising critical 
point to be at J/T c = 0.2216544(6) (i.e. T c /J « 4.51), and of course we already know that 
H c = 0. What we now want to figure out is whether or not we can infer the location of this 
critical point by carrying out the expansion described above. 



III. RESULTS 

Monte Carlo simulations were performed on 3D simple cubic lattices using the Cluster 
Algorithm [2H]- As mentioned in the introduction, we are restricted on the phase diagram 
to lines of constant background magnetic field strength, denoted by H Q . We specifically 
choose to work along the lines of H = 0.1J and H = 0.01J, since at these distances we 
are able to obtain good statistics, while remaining in the critical region (when near T c ). At 
distances too close to the critical point the correlation length grows rapidly, hence much 
larger volumes would require consideration. Whereas, if one is too far from the critical point 
(Ho = J, say), critical behavior is not observed. For the values of H that we considered, it 
was sufficient to consider sizes in the range 4 < N < 24. 
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Now that we have chosen the values of Ho around which we will base our analysis, we need 
to determine the minimal lattice sizes (for each of our the two different values of H Q ) which 
sufficiently converge to the thermodynamic limit. Due to the divergence of the correlation 
length at the critical point, one would expect that obtaining reasonable data for H = 0.01J 
requires a larger (perhaps significantly) lattice than for H = 0.1J. 

Such convergence is best observed qualitatively by plotting Ki^. CIS db function of T over a 
range of volumes (shown in Fig. 3(a) for Ho = 0.1J and Fig. 3(b) for Ho = 0.01J). These 
graphs suggest that it is sufficient to work with N = 8 and N = 20 for Ho = 0.1J and 
Ho = 0.01J respectively, based on qualitative observation that the thermodynamic "limit" 
is achieved whenever any two subsequent graphs overlap. 

More to the point, the thermodynamic limit is effectively reached when the volume de- 
pendence of higher order cumulants becomes washed away in statistical uncertainty. Said 
differently, it becomes impossible to accurately distinguish volume dependence from noise. 
This effect is shown in Figs. 4(a) and 4(b) for Kg, where it is made very clear for Ho = 0.01J. 
The increasing uncertainty eventually renders the analysis inconclusive at sufficiently large 
volumes. 

Based on Figs. |3(a)| and 4(a)| it may seem a bit overzealous to regard N = 8 as the 
thermodynamic limit. However, at the expense of relatively small finite-size volume effects, 
it turns out to be very statistically advantageous to work with N = 8 over say N = 9; 
this is further discussed in the next section. It should be noted that measurements of the 
correlation length along our values of H (see Appendix A) further validate the choices of 
N = 8 and N = 20, since it turns out that these lattices are many correlation lengths across. 

Henceforth, if we are discussing H = 0.1J results it should be assumed that the data 
originated from an N = 8 lattice and likewise N = 20 for Hq = 0.01J, unless it is stated 
otherwise. Moreover, it should be assumed that any measured quantity is obtained from 
averaging over Af = 10 7 uncorrected samples, with estimates of uncertainty calculated using 
the Jackknife Method [30]. 

The results of our analysis are plotted over Figs. [5] to 11 Fig. [5] illustrates the crossover 
from the magnetized to the unordered phase (ki is simply the magnetization per unit vol- 
ume). The Taylor coefficients (higher order k^s, normalized to unity) are shown in Fig. 
|6j The Taylor series constructed out of these coefficients are depicted in Figs. 7] and [TT 
Finally, r n (from which we extrapolate radii of convergence) are shown in Figs, p to 10 



IV. ANALYSIS AND DISCUSSION 

Now that the data has been presented, the first thing that we should try to do is infer 
(guess) the location of T c . Subsequently, we will then calculate H c and specify the location 
of the critical point (assuming that everything goes according to plan). 

Monte Carlo simulations of the Ising Model are drastically simpler than those of QCD. 
That being said, the "restricted" extent of our data set is not due to any sort of numerical 
limitations. If we truly desired we could easily a) obtain data from larger sample sizes, 
b) sample the temperature axis more frequently and c) obtain even higher order Taylor 
coefficients. Rather, as discussed in the introduction, we are on purpose trying to somewhat 
mimic the limitations facing any realistic study of QCD. 

In what follows, we will try to apply to same logic in locating the Ising critical point as 
that which has been used for QCD. Therefore, a bit of a synopsis of the analysis in [21] 
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Ql p q2 p 

would be in order. The authors focus mainly on the Taylor expansions of » Q 2 CD and 75 — ¥^ 
up to sixth order 2 . By definition, expanding in fi = /xb/3 and \ij = 0, 

^ 2 — N ,0 1 2 2 1 4 4 1 6 



2 P QCD (T,/x) = + X> 2 + X> 4 + X> b + (4-1) 

where x°b = X2,o and, for instance x% = (X6,o+4x5,i+'' 7 X4,2+4x3 i 3)/4! with similar expressions 
for the other Xs' s - The critical region is highlighted by the formation of sharp peaks in the 
Xb's, which are observed in the vicinity of T c / Q . Furthermore, at T < 0.95T c / o , x% changes 

q2 p 

sign. By plotting the Taylor series of g^2 CD ; criticality is most evident at T = 0.95T c / o , 
since this temperature corresponds to the most rapidly diverging graph. Subsequently, the 
distance to the CEP if inferred from an extrapolation of the series' radius of convergence. 
Hence, in our case, the analysis will proceed in the following order, 

• locate the critical region 

• observe the range over which the Ising coefficients exhibit the correct alternating sign 
behavior indicative of a divergence along the negative H axis 

• estimate the exact location of T c by plotting x's Taylor series 

• finally, determine the location of the critical point from our extrapolations of r n . 

Since we are not on the H = axis, nor are we in the infinite volume limit, the curves 
in Fig. [6] depict a crossover rather than a phase transition. By observing the range over 
which k 8 exhibits non-trivial behavior, we can safely say that T c j J is bracketed between 
[4.25, 5.00] for H = 0.1J and [4.45,4.65] for H = 0.01J. 

Continuing on, the singular point that we are looking for lives on the negative H axis of 
our Taylor series. In other words, it is a requirement for our Taylor coefficients to alternate 
in sign order by order. Looking back at Fig. |6j we see that at T = 4.80 J and T = 4.60 J for 
Hq = 0.1J and Ho = 0.01J respectively the 7 th order coefficient changes sign, thus breaking 
the pattern. From this observation we can refine our upper bounds on T c . 

At the order that we are working with, there is a whole range (rather than a single point) 
in T where the coefficients are observed to alternate in sign. Hence, this property alone 
is not sufficient to specify T c . If we were to plot the full Taylor series for the magnetic 
susceptibility x, we would observe a genuine singularity at the critical point T c (assuming 
that we are in the thermodynamic limit). Aside from being at finite volume, we do not have 
the full series, so the best that we can do is plot the partial sum of the first 6 terms (we only 
have 6 rather than 8 terms since we wish to plot x rather than the free energy F, analogous 
to the situation in QCD). 

The rate at which this polynomial diverges is directly attributable to the values of the 
Taylor coefficients, i.e. the temperature at which this polynomial diverges most rapidly 
should give a somewhat clearer indication of the location of T c . The partial sum of x' s 
Taylor series up to order n is given by 



i=i 

2 Note that even when m = m u = m c i and [i = fi u = [id, X2.0 does not necessarily equal xi,i (flavour 
symmetry only requires X2,o = Xo,2)- Simply, xi,i ~ Det'Det' while X2.0 ~ Det"Det. 
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plotted in Fig. [Tj Uncertainties, when depicted, are due to the standard error 
of the highest order Taylor coefficient; in this case that would be Kg. When Ho = 0.1J, 
the T = 4.50J, T = 4.55J and T = 4.60J graphs, aside from being indistinguishable, 
diverge most rapidly. From the plots we expect that 4.50 < T c /J < 4.60. Furthermore, the 
H = 0.01J plot suggests that T c /J ~ 4.50. 

Now that we are led to believe that T c j J = 4.55 ± 0.05 (or least that T c is in the vicinity 
of 4. 55 J), we will try to locate H c . In the infinite volume limit, 

H c = H -T croo (T c ), (4.3) 

where r^T) is the radius of convergence as defined earlier, which in general is temperature 
dependent, and most importantly r^iTc) is finite due to non-analyticity at the critical point. 
However, for any finite volume, the partition function is analytic everywhere on the phase 
diagram, hence the radii of convergence that we measure should always diverge when n — > oo. 
So, in fact, when we say that we are extrapolating r n to (i.e. taking the limit defined by 



Eq. 2.14), what we really mean is that we expect to observe a range in values of n where ry, 



forms a plateau. The value of r n along this plateau is used to determine H c via Eq. 4.3 



In Fig. [8j we observe exactly this sort of behavior. r n (T) forms a plateau for some 
(but not all) values of T in the vicinity of T c (shown in greater detail in Fig. 10). With 
Ho = 0.01J, a plateau forms when T = 4.50J, and it continues all the way up to n = 7. 
When we go to T = 4.55J, we see that by n = 7 the plateau has started to tilt upwards. 
The H = 0.1J graphs exhibit similar behavior, see Fig. [9j 

For both values of H , r n has the smallest error bars when T = 4.50 J; furthermore, the 
plateaus that form at this temperature extend all the way up to n = 7. Hence, our most 
reliable estimate of H c can be obtained from the T = 4.50J graphs. Within error bars, 
-0.005 < HJJ < 0.015 from H = 0.1J and -0.002 < HJJ < 0.002 from H = 0.01J. 

Thus, from the analysis centered about H = 0.1J, we can conclude that 4.50 < TJ J < 
4.60 and —0.005 < H c /J < 0.015, a finding which is in perfect agreement with the known 
location of the Ising critical point. Therefore, we have observed that this method is able to 
effectively locate the Ising Model critical point. 

In working with the Ising model, we were not limited numerically, hence we carried out 
this analysis to 8 th order in the cumulant expansion without much difficulty. In doing so, we 
were able to make a few observations that may be of relevance to a study revolving around 
QCD. 

Foremost, our best data came from the N = 8 and iV = 20 lattices, which were the 
minimum sizes that could be considered "thermodynamic." As we discussed, the subtraction 
of disconnected pieces from high order cumulants leads to a very rapid deterioration with 
volume in the statistics. Therefore one must make a careful balancing act between large 
volumes to control finite-volume effects, and the loss of statistics in high-order cumulants 
with large volume. 

It is also interesting to see how the effectiveness of the method changes as one is limited to 



lower orders in the Taylor series. Plots of x^ an d X^ are shown in Fig. 11 When limited 
to 6 th order, it would appear that 4.60 < TJ J < 4.65. This is a slight overestimate, but the 
real issue here is that it becomes impossible to determine H c . In Figs. |8j [9] and 10 we see 
that r n only forms a plateau for n > 5. At 8 th order, we extrapolate r n solely from rs, r& and 
r-j. However, when restricted to 6 th order, no such extrapolation is possible, since you would 
only know r^. Recall that, since the expansion in QCD in terms of the chemical potential 
only contains even powers in fi, the 4'th order cumulant in the Ising model is equivalent to 
the n 8 term in QCD. 
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V. CONCLUSION 



We have tested the prospects for the Taylor Series Technique for locating a critical end- 
point in a phase space which cannot be explored directly. By finding the cumulants in terms 
of the unexploreable direction, one can extrapolate the radius of convergence and type of 
convergence to find a critical temperature and other parameter (chemical potential in QCD). 
In the "toy" case where we confine ourselves to the Ho = 0.1J line of the 3D Ising model 
phase diagram, we calculated the Taylor coefficients in an expansion of the free energy in 
H. By analyzing the properties of these Taylor coefficients, as well as extrapolating the 
radius of convergence of the series which they form, we were led to conclude that the critical 
point lives in the range 4.50 < TJJ < 4.60 and -0.005 < HJ.J < 0.015. This finding is 
in agreement with the known result that T c /J pa 4.51J and H c /J = 0. Therefore, we find 
that the Taylor series method can be applied successfully to find a critical endpoint outside 
of the range where direct lattice studies are conducted. In this sense, the findings of this 
analysis offer a positive outlook for the applicability of the method to QCD. 

However, we found that a successful extrapolation of the critical point requires high orders 
in the cumulant expansion, as well as excellent control of errors. In the case we considered, 
an accurate determination of the distance to the critical endpoint required the use of 8 terms 
in the cumulant expansion; using only 4 terms (the number currently available in the case 
of QCD) is not enough to determine the distance to the critical point accurately. 



Appendix A: Correlation Lengths 



In Section III we concluded that N = 8 and N = 20 lattices sufficiently approximate 



the thermodynamic limit when H = 0.1J and H = 0.01J respectively. This was based 



primarily on observations of the volume dependence of Kg (see Figs. 4(a) and 4(b) ). 

Alternatively, we could have reached this conclusion by calculating correlation lengths, 
denoted by £. For our purposes it was sufficient to qualitatively observe the convergence of 
Kg, however, it is still of interest to examine the thermal dependence of £ for both background 
field strengths. This will serve as a non-trivial check that it was sensible to base our findings 
on N = 8 and N = 20 lattices, i.e. we should observe that iV 3> £. 

£ can be calculated by measuring the exponential fall-off of the two-point function G(r) = 
(a(x)a(O)). The tail of G(r) was subsequently fitted to 

r 

G(r) = + M 2 (Al) 

since G exhibits the correct asymptotic behaviour, rj was fixed at 0.025, in accordance with 
the results found in 



Based on the results of this fit, £ as a function of temperature is plotted in Fig. [12 
These graphs are somewhat noisy since measurements of G(r) were obtained from M = 10' 
uncorrelated lattices. However, what we observe is that £ ~ 1 when Hq = 0.1J, while £ 
jumps to near 4 when H = 0.01J, justifying the necessity to go to a larger volume. 
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FIG. 3: Volume dependence of K4 indicating convergence to the thermodynamic limit. 
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FIG. 4: Volume dependence of Kg at T = 4.60J (note the volume normalization, the quantity 
plotted is actually V 1 Kg). 
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(a) H = 0.1J (b) H = 0.01J 

FIG. 5: Temperature dependence of the magnetization per unit volume, showing the crossover 
between the ordered / unordered phases. 
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FIG. 6: Temperature dependence of the higher order Taylor coefficients, normalized so that multiple 
curves can be depicted on each plot. Even and odd coefficients are separated for clarity. 
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FIG. 7: as defined by Eq. 4.2 The leftmost graphs depict the progression of throughout 
the critical region (uncertainties are not depicted). The rightmost plots contain the most rapidly 
divergent graphs indicating that one is very near T c . Upper and lower curves correspond to uncer- 
tainty in Kg. When Hq = 0.1J the T = 4.50J, T = 4.55J and T = 4.60J are indistinguishable. 
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FIG. 9: Close-up of the H = 0.1J graphs in Fig. ffl 
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FIG. 10: Close-up of the H = 0.01J graphs in Fig. M 
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FIG. 11: an d X*' 5 ''- When restricted to lower orders, we see that one would be inclined to 
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FIG. 12: Correlation length £, obtained by fitting Eq. Al to (cr(r)cr(O)), illustrating why it was 
necessary to go to N = 20 for Hq = 0.01J. 
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